В данной работе изпользуются данные о ценности жилья в городе Бостон в 1970-1980-х годах (пакет MASS).
Основной акцент будет сделан на то, как средняя стоимость домов (переменная medv, измеренная в тысячах долларов), занимаемых владельцами, зависит от различных факторов.
Авторы работы измерили значения 14 различных параметров (включая medv) для 506 домов.
Цель — построить линейную модель для предсказания стоимости домов.
Используемые обозначения:
Примечание в конце анализа я осознал, что переменная rad, в целом, является фактором, но у него большое количество уровней, поэтому я решил не переделывать, подумав, что это не сильно повлияет на анализ.:)
Давайте посмотрим, как выглядят данные.
house_cost <- Boston
house_cost$chas <- factor(house_cost$chas)
house_cost$chas <- relevel(house_cost$chas, ref = '0')
str(house_cost) # Все ли нормально открылось
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
plot(house_cost) # скатерплот чтобы посмотреть структуру данных
Для начала проведем стандартизацию переменных. Также на первом этапе не будем учитывать взаимодействие предикторов.
house_cost_scale <- as.data.frame(sapply(house_cost[-4], scale))
house_cost_scale$chas <- house_cost$chas
house_cost_scale$chas <- relevel(house_cost_scale$chas, ref = '0')
Теперь построим полную линейную модель
mod_full <- lm(medv ~ ., data = house_cost_scale)
summary(mod_full)
##
## Call:
## lm(formula = medv ~ ., data = house_cost_scale)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.69559 -0.29680 -0.05633 0.19322 2.84864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.020206 0.023835 -0.848 0.396976
## crim -0.101017 0.030737 -3.287 0.001087 **
## zn 0.117715 0.034811 3.382 0.000778 ***
## indus 0.015335 0.045871 0.334 0.738288
## nox -0.223848 0.048126 -4.651 4.25e-06 ***
## rm 0.291056 0.031928 9.116 < 2e-16 ***
## age 0.002119 0.040430 0.052 0.958229
## dis -0.337836 0.045666 -7.398 6.01e-13 ***
## rad 0.289749 0.062813 4.613 5.07e-06 ***
## tax -0.226032 0.068912 -3.280 0.001112 **
## ptratio -0.224271 0.030796 -7.283 1.31e-12 ***
## black 0.092432 0.026662 3.467 0.000573 ***
## lstat -0.407447 0.039378 -10.347 < 2e-16 ***
## chas1 0.292128 0.093679 3.118 0.001925 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.516 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
Ни одно значение не превышает условного порога в 2 единицы. Влиятельных наблюдений нет (Но, если честно, то видно, что есть ряд наблюдений, которые могут смещать оценки коэффициентов модели, так как для них значения расстояния Кукa больше).
ggplot(mod_full_diag, aes(x = 1:nrow(mod_full_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red") +
xlab('Номер наблюдения') +
ylab('Расстояние Кука')
График распределения остатков выглядит не очень хорошо. Достаточно большое количество наблюдений находятся за пределами двух стандартных отклонений, а также отчетливо виден паттерн в остатках. Все это свидетельствует о наличии нелинейности во взаимосвязи, а также о гетероскедастичности (непостоянство дисперсии).
gg_resid <- ggplot(data = mod_full_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_hline(yintercept = 0) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red") +
xlab('Предсказанное значение') +
ylab('Oст. модели')
gg_resid
Квантильный график выглядит не очень хорошо, сказать, что стандартизованные остатки распределены нормально, нельзя.
qqPlot(mod_full_diag$.stdresid, xlab = 'Квантили нормального распределения', ylab = 'Квантили распределения остатков модели')
## [1] 369 372
Если обратиться к структуре датасета, то можно заметить наличие корреляций между предикторами. Поэтому в данном случае проверка на мулькиколлинеарность особенно актуальна.
Один из способов проверки модели на мультиколлинеарность предикторов это использование VIF (variance inflation factor).
Если предиктор имеет значение VIF выше 2, то его следует исключать из модели.
Если таких предикторов несколько, то прменяется пошаговый алгоритм: расчет VIF, удаление предиктора с максимальным VIF, расчет VIF для обновленной модели до тех пор, пока все значения не будут меньше порога.
vif(mod_full)
## crim zn indus nox rm age dis rad
## 1.792192 2.298758 3.991596 4.393720 1.933744 3.100826 3.955945 7.484496
## tax ptratio black lstat chas
## 9.008554 1.799084 1.348521 2.941491 1.073995
В даном случае легко можем видеть наличие мультиколлинеарности в нашей модели. В дальнейшем ряд предикторов надо будет убрать из модели.
Построим график предсказаний стоимости от переменной, которая обладает наибольшим по модулю коэффициентом. В данном случае это переменная lstat.
Мы не можем одновременно учесть изменчивость всех предикторов. Поэтому как правило выбирается один, который интересует вас больше всего (в нашем случае lstat). И относительно него создается тестовый датасет, где целевой предиктор принимает значения от минимального до максимального, а все остальные предикторы представлены своими средними значениями.
test_data <- data.frame(
lstat <- seq(min(house_cost_scale$lstat), max(house_cost_scale$lstat), length.out = 90),
crim <- rep(mean(house_cost_scale$crim), 90),
zn <- rep(mean(house_cost_scale$zn), 90),
indus <- rep(mean(house_cost_scale$indus), 90),
nox <- rep(mean(house_cost_scale$nox), 90),
rm <- rep(mean(house_cost_scale$rm), 90),
age <- rep(mean(house_cost_scale$age), 90),
dis <- rep(mean(house_cost_scale$dis), 90),
rad <- rep(mean(house_cost_scale$rad), 90),
tax <- rep(mean(house_cost_scale$tax), 90),
ptratio <- rep(mean(house_cost_scale$ptratio), 90),
black <- rep(mean(house_cost_scale$black), 90),
chas <- rep('0', 90)
)
Predictions <- predict(mod_full, newdata = test_data, interval = 'confidence')
MyData <- data.frame(test_data, Predictions)
Pl_predict <- ggplot(MyData, aes(x = lstat, y = fit)) +
geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
geom_line() +
ggtitle("Множественная модель") +
xlab('Низкий статус населения, sd') +
ylab('Предсказанное значение цены за дом, sd')
Pl_predict
Полученная модель далека от идеальной, как было показано выше, также мы не убрали предикторы, которые незначимо вдияют на зависимую переменную. Она уже может делать, какие-то предсказания, но ниже мы попытаемся сделать ее лучше.
Вернемся к нашей полной модели и вспомним результаты проверки на коллинеарность.
vif(mod_full)
## crim zn indus nox rm age dis rad
## 1.792192 2.298758 3.991596 4.393720 1.933744 3.100826 3.955945 7.484496
## tax ptratio black lstat chas
## 9.008554 1.799084 1.348521 2.941491 1.073995
Будем пошагово удалять из модели по одному предиктору с наибольшим vif до тех пор, пока все значения не будут меньше порога (2).
mod2 <- update(mod_full, .~. - tax)
vif(mod2)
## crim zn indus nox rm age dis rad
## 1.791940 2.184240 3.226015 4.369271 1.923075 3.098044 3.954446 2.837494
## ptratio black lstat chas
## 1.788839 1.347564 2.940800 1.058220
mod3 <- update(mod2, .~. - nox)
vif(mod3)
## crim zn indus rm age dis rad ptratio
## 1.785343 2.183394 2.872809 1.904013 2.875130 3.641492 2.533616 1.598944
## black lstat chas
## 1.339554 2.927273 1.057571
mod4 <- update(mod3, .~. - dis)
vif(mod4)
## crim zn indus rm age rad ptratio black
## 1.765881 1.758636 2.517520 1.879925 2.423551 2.507024 1.530992 1.339553
## lstat chas
## 2.926111 1.056840
mod5 <- update(mod4, .~. - lstat)
vif(mod5)
## crim zn indus rm age rad ptratio black
## 1.727110 1.751140 2.489395 1.307312 2.038976 2.497242 1.529761 1.304899
## chas
## 1.053402
mod6 <- update(mod5, .~. - nox)
vif(mod6)
## crim zn indus rm age rad ptratio black
## 1.727110 1.751140 2.489395 1.307312 2.038976 2.497242 1.529761 1.304899
## chas
## 1.053402
mod7 <- update(mod6, .~. - rad)
vif(mod7)
## crim zn indus rm age ptratio black chas
## 1.372553 1.738620 2.229252 1.281667 2.026045 1.369726 1.251785 1.052519
mod_good <- update(mod7, .~. - indus)
vif(mod_good)
## crim zn rm age ptratio black chas
## 1.342583 1.680416 1.216701 1.668190 1.345528 1.209887 1.043732
Итоговая модель этого шага (в нестандартизованном виде):
medv = -7.06 - 0.11 * crim - 0.003 * zn + 7.06 * rm - 0.04 * age - 0.93 * ptratio + 0.015 * black
medv = -3.54 - 0.11 * crim - 0.003 * zn + 7.06 * rm - 0.04 * age - 0.93 * ptratio + 0.015 * black
Мы можем оставить модель в таком виде, а можем попробовать оставить только те предикторы, которые значимо влияют на стоимость домов (большое число предикторов обычно приводит к оверфитингу модели).
В этой работе я буду использовать алгоритм отбора предикторов backward selection (он же backward elimination). В качестве критерия отбора буду использовать частный F-тест.
drop1(mod_good_unstd, test = 'F')
## Single term deletions
##
## Model:
## medv ~ crim + zn + chas + rm + age + ptratio + black
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 15509 1747.9
## crim 1 323.2 15832 1756.3 10.3772 0.0013592 **
## zn 1 1.6 15511 1745.9 0.0506 0.8220598
## chas 1 386.8 15896 1758.3 12.4206 0.0004638 ***
## rm 1 10225.3 25735 2002.1 328.3350 < 2.2e-16 ***
## age 1 411.9 15921 1759.1 13.2258 0.0003049 ***
## ptratio 1 1514.7 17024 1793.0 48.6365 9.843e-12 ***
## black 1 781.1 16290 1770.7 25.0805 7.646e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod_good_unstd <- update(mod_good_unstd, .~. - zn)
drop1(mod_good_unstd, test = 'F')
## Single term deletions
##
## Model:
## medv ~ crim + chas + rm + age + ptratio + black
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 15511 1745.9
## crim 1 329.2 15840 1754.5 10.591 0.0012133 **
## chas 1 390.0 15901 1756.5 12.546 0.0004342 ***
## rm 1 10413.2 25924 2003.8 335.006 < 2.2e-16 ***
## age 1 514.5 16025 1760.4 16.553 5.498e-05 ***
## ptratio 1 1605.3 17116 1793.7 51.644 2.440e-12 ***
## black 1 780.3 16291 1768.8 25.103 7.559e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Таким образом, мы оставили только предикторы со значимым влиянием на зависимую переменную.
Мы уже смотрели на ряд графиков для диагностики нашей модели, но помимо них для множественных моделей необходимо строить графики от предикторов, не вошедших в модель. В нашем случае видно, что в модели, в общем-то, нет неучтенных зависимостей, соответственно, можно не возвращать предикторы, убранные из модели на предыдущих этапах анализа
res1 <- gg_resid + aes(x = zn)
res2 <- gg_resid + aes(x = indus)
res3 <- gg_resid + aes(x = nox)
res4 <- gg_resid + aes(x = dis)
res5 <- gg_resid + aes(x = rad)
res6 <- gg_resid + aes(x = tax)
res7 <- gg_resid + aes(x = lstat)
grid.arrange(res1, res2, res3, res4, res5, res6, res7, nrow = 4)
mod_good_unstd_diag <- data.frame(fortify(mod_good_unstd), house_cost[, c(1, 4, 6, 7, 11, 12)])
ggplot(mod_good_unstd_diag, aes(x = 1:nrow(mod_good_unstd_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red") +
xlab('Номер наблюдения') +
ylab('Расстояние Кука')
gg_resid <- ggplot(data = mod_good_unstd_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_smooth(method = "lm") +
geom_hline(yintercept = 0) +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red") +
xlab('Предсказанное значение') +
ylab('Стандартизованные остатки модели')
gg_resid
qqPlot(mod_good_unstd_diag$.stdresid, xlab = 'Квантили нормального распределения', ylab = 'Квантили распределения остатков модели')
## [1] 369 372
Анализ показал, что, если неприятности и исправились, то не очень сильно. Модель вышла не идеальная. Может быть, в данном случае помогла бы какая-то трансформация переменных (зачастую не желательна, так как дальше сложно интерпретировать результаты) или применение других статистических методов, о которых, вероятно, речь пойдет в будущих проектах. Также попытки учесть какие-то взаимодействия переменных внутри итоговой модели не привели к успеху, наоборот график расстояний Кука становился только хуже (не привожу здесь). В первом приближении мы можем попытаться предсказать, какие аспекты надо улучшить, чтобы максимизировать цену за дом.
Итоговая модель:
medv = -7.06 - 0.11 * crim + 7.05 * rm - 0.04 * age - 0.92 * ptratio + 0.015 * black
medv = -3.54 - 0.11 * crim + 7.05 * rm - 0.04 * age - 0.92 * ptratio + 0.015 * black
Мы можем видеть 6 параметров, которой так или иначе влияют на цену домов. Вот несколько замечаний, которые, на мой взгляд, должны максимизировать цену за продаваемый дом:
Таким образом в идеальном районе:
Примерная цена за дом при выполнении вышеуказанных условий может составить 50.1 тысяч долларов.